 function [ysim,stmat,shocks]=simgensys_data(G1,C,impact,SDX,zmat,tvec,Y0)
%function [ysim,stmat]=simgensys(shocks,G1,C,impact,zmat,Y0);
% Simulate the path of the variables given the output of gensys.m  
% Inputs 
% ------
% shocks    [NX x T ] vector of shoclk
% G1        Output from gensys 
% C         Constant, Output from gensys
% impact    Output from gensys
% zmat      (optional)  Matrix multiplies the state vector resulting in a matrix of 
%           observables 
%           If not defined, set to the identity matrix 
% Y0        (optional) first state 
%           If not defined, set to zero
%
% Alejandro Justiniano 5/27/03 
% Modified 5/20/04 Added the constant
% Oct 2 2007  Matches what would obtain from KFILTER if fed the 
% simulated shocks here 
% ------------------------------------------------------
if ~isequal(size(C)',[size(G1,1);1] ) ; 
    error('Constant AR matrix dimension mistmatch') 
end 
if length(tvec)~=2;error('TVEC must be 2x1');end
dispaj('Discard first ',tvec(1),' keep ',tvec(2));
nx=size(SDX,1); 
ndr=sum(tvec); 
shocks=(randn(ndr,nx)*(SDX'))'; 
ysim=zeros(size(zmat,1),ndr); 
stmat=zeros(size(G1,1),ndr); 
if nargin < 7; 
    stmat(:,1)=zeros(size(G1,1),1); 
else 
    stmat(:,1)=Y0;
end
ysim(:,1)=zmat*stmat(:,1); 
%disp(['Begin simulation for ',num2str(T),' periods '] ); 
for ii=2:ndr;     
    stmat(:,ii)=G1*stmat(:,ii-1)+impact*shocks(:,ii);
    ysim(:,ii)=zmat*(stmat(:,ii)+C);  
end 
disp('Truncate'); 
ysim=ysim(:,tvec(1)+1:end)'; 
stmat=stmat(:,tvec(1)+1:end)'; 
shocks=shocks(:,tvec(1)+1:end)'; 